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ABSTRACT 

We introduce a new project to understand helium reionization using fully coupled A^-body, hydrody¬ 
namics, and radiative transfer simulations. This project aims to capture correctly the thermal history 
of the intergalactic medium (IGM) as a result of reionization and make predictions about the Lyman-a 
forest and baryon temperature-density relation. The dominant sources of radiation for this transition 
are quasars, so modeling the source population accurately is very important for making reliable pre¬ 
dictions. In this first paper, we present a new method for populating dark matter halos with quasars. 
Our set of quasar models includes two different light curves, a lightbulb (simple on/off) and sym¬ 
metric exponential model, and luminosity-dependent quasar lifetimes. Our method self-consistently 
reproduces an input quasar luminosity function given a halo catalog from an A^-body simulation, and 
propagates quasars through the merger history of halo hosts. After calibrating quasar clustering using 
measurements from the Baryon Oscillation Spectroscopic Survey, we find that the characteristic mass 
of quasar hosts is Mh ~ 2.5 x 10^^ h~^MQ for the lightbulb model, and Mh ~ 2.3 x 10^^ h~^MQ 
for the exponential model. In the latter model, the peak quasar luminosity for a given halo mass is 
larger than that in the former model, typically by a factor of 1.5-2. The effective lifetime for quasars 
in the lightbulb model is 59 Myr, and in the exponential case, the effective time constant is about 15 
Myr. We include semi-analytic calculations of helium reionization, and discuss how to include these 
quasars as sources of ionizing radiation for full hydrodynamics with radiative transfer simulations in 
order to study helium reionization. 

Keywords: cosmology: theory — intergalactic medium — large-scale structure of the universe — 
quasars: general 


1. INTRODUCTION 

Helium reionization is an important epoch in the Uni¬ 
verse’s history, and the most recent large-scale transition 
of the intergalactic medium (IGM). During the epoch of 
hydrogen reionization, the first stars and galaxies emit¬ 
ted photons capable of ionizing hydrogen and singly ion¬ 
izing helium (whose ionization energies are 13.6 and 24.6 
eV, respectively). However, the spectra of these first 
sources did not contain a sufficient number of high-energy 
photons capable of doubly ionizing helium, which re¬ 
quires a much larger ionization energy (54.4 eV). Conse¬ 
quently, helium was predominantly singly ionized follow¬ 
ing hydrogen reionization until a burst of quasar activity 
at redshifts 6 > 0 > 2. Quasars are thought to be the 
first objects to emit an appreciable number of photons 
capable of doubly ionizing helium. However, because the 
birth of quasars requires additional time for structure to 
form and sufficient mass to assemble inside dark matter 
halos, this period of evolution occurs later in the Uni¬ 
verse’s history. 


Recent and upcoming efforts to look for quasars include 
the Baryo n Oscillation Spectr oscopic Survey (BOSS) of 


SDSS-Hl ( 
of the Sul 

Dawson et al. 12013 1, the Hyper Suprime Cam 

oaru telescope (Kashikawa et al. 2015), and 

DESl (ISchlegel et al. 20111. I’here are currently abont 

420,000 ui 
number pi 

lique quasar objects (Fiesch |2015 1, with this 

ejected to increase by an order of magnitude 


after the conclusion of the next generation of experi¬ 
ments. This rich set of observations allows us to char- 
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acterize quasars to an unprecedented level of accuracy, 
and better characterize their properties. This is espe¬ 
cially true at high redshift (^ ^ 6), where there are cur¬ 
rently few observations. Determining quasar properties 
at high redshifts is helpful for understanding the growth 
of structure, as well as providing observations of reion¬ 
ization through measuring their absorption spectra. 
Observations have shown that quasar activity peaks 


betwe en 2 < x < 3 (iWarren et al.||1 994[ Schmidt et ah, 
1995 1. The Gunn-Peterson trough ([Cunn &: Reterson 


1965 1 of he l ium has been dete c ted at z > 3 (|jakobsen 
et al.||1994t |Zheng et al.||2008| |Syphers & Shull||2U14p, 
implying that some fraction of helium was still present 
as He II at these redshifts. Helium absorption then tran¬ 
sitions to becoming patchy, with extended regions of ab- 
sorption and transmis sion in the He ii Lyman-a forest 
dReimer s et al.| [19971), and seems to be completed by 
z ^ 2.7 ( Dixon fc Eurlanetto||2009 Worseck et al.|20jT ), 
which coincides with the peak in quasar activity. How- 
ever, to observe the Gunn-Peterson trough of He ii, the 
sight line must be free of any intervening Lyman-limit 
systems. This means that the number of observations 
for these measurements is rather small (of 0(10)). 

When discussing helium reionization, it is important 
to understand the properties of the ionization sources, 
such as quasars’ lifetimes and light curves. On the the¬ 
oretical side of the problem, there are some predictions 
for quasar properties, but also a fair degree of uncer¬ 
tainty. By treating quasars as accretion disks around 
super-massive black holes (SMBHs), one can show that 
the maximal conversion efficiency e for converting mass 
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to luminosity is e ^ 0.3 (Thorne 1974). Further, for 
matter accreting onto an SMBH at the Eddington limit 
(Eddington 1926), one obtains an exponential increase 
in mass and luminosity with a characteristic time scale 
(called the Salpeter e-folding time) of r = 45 Myr for 
e = 0.1 (|Salpeterjl964||Wyithe k. Loeb|2003|). Cosmolog¬ 
ical simulations that seek to capture the relationship be¬ 
tween quasars and their galaxy hosts have treated quasar 
activity as being th e result of a major-merger event be¬ 
tween two galaxies (Springel et al. 2005 Hopkins et al. 


2006 2008), or a cold-how accretion ot gas onto the cen¬ 
tral SMBH (Di Matteo et al. 2012). However, there is 
no definitive evidence that quasars accrete exclusively at 
the Eddington limit, or are limited to a single episode of 
highly luminous activity. 

Observations can also help us understand the physics 
of quasars, though typically at larger scales than theory 
or simulation. Since the entire rise and fall of quasar 
number density spans a time of roughly 1 0^ years, the 


quasar lifetime must be shorter than this (Osmer 2004 
p. 324). At the other extreme, observations of the quasar 
proximity zone show that qua sar lifetimes should be at 
least 10® years (Martini 2004 p. 169). This time scale 


corresponds to the photoionization timescale of relatively 
high-density neutral hydrogen systems observed to be 
ionized in the IGM, and so the lifetime of the quasar 
must be at least this long in order to maintain the highly 
ionized level of these systems observed in the Lyman-a 
forest. Eurther constraints are difficult to obtain, and 


ing measurements (e.p., Porciani et al.||2004 

Porciani & 

Norberg 

2U06 White et al. 

2U12|). Estimates made us- 


that are 10-100 Myr, with most values being ~30 Myr, 
which is comparable to the Salpeter e-folding time. Fur¬ 
ther, there are few definitive constraints on qu asar light 
curves (though see Hopkins fc Hernquist||2009 ). 


For the universal populations ot quasars, the major 
pieces of data are their number density as a function of 


tion (QLF) (t>{L.z). e.c 

., Schmidt & Green 1983 Boyle 

et al.| 20001 Ross et al. 

2U13D, and their spatial cluster- 

ing (KJutram et al.|2U03 

Porciani et al.|2004 

White et al. 

2U12p. These observations can constrain scaling relations 


2013) , or used to calibra te subgrid models for simulations 


{e.g., Feng et al. 2014). However, as mentioned above, 
the properties ot individual quasars are difficult to ex¬ 
tract from these observations, due to degeneracies. The 
imposed constraints are typically weak, and only provide 
order-of-magnitude precision. 

Cosmological simulations are an ideal tool for further¬ 
ing our knowledge about this portion of the universe’s 
history. Helium reionization leaves a lasting impression 
on the thermal history of the IGM: the relative hardness 
of quasar spectra means that there is a large degree of 
photoheating of the IGM while reionization is occurring. 
Thus, it is important to include hydrodynamics in simu¬ 
lations, in order to include the effects of baryonic physics. 
Additionally, due to the relatively long mean free path of 
far-UV and soft X-ray photons when looking at helium 
reionization, it becomes important to include radiative 
transfer calculations in simulations. Thus, semi-analytic 
calculations that assume a sharp reionization front are 


typically poor approximations of the physical situation. 
Even ID radiative transfer codes are not realistic enough 
to calculate the inhomogeneous reionization process, es¬ 
pecially when reionized regions begin to overlap. Due to 
the highly biased nature of quasar sources, this is typ¬ 
ically early in the reionization process. Therefore, 3D 
radiative transfer calculations are essential for capturing 
the complicated physics of helium reionization. As men¬ 
tioned earlier, the large degree of thermal heating argues 
for simulations in which the hydrodynamics calculations 
are coupled to the radiative transfer ones. This work 
builds on and extends previous investigations o f helium 


reionization, wh i ch either were sem i-numerical (Furlan- 
etto fc Oh||2008 Dixon et alj|2014) or applied radiative 


transfer m post-processing ( McQuinn et al. 2009, 2011 
Compostella et ^|2013 2015(7” 


(Jur approach to helium reionization uses simulations, 
with Af-body, hydrodynamics, and radiative transfer 
solved simultaneously. An essential first step of this cal¬ 
culation is to understand the sources of reionization, and 
ensure that their properties match the observations as 
nearly as possible. To this end, we use the observed QLF 
from the SDSS and the COSMOS survey across various 
redshift epo chs (Masters et al.|2012 McGreer et al.|2013 
Ross et ^|2013 hereafter lVfl2, lVfl3, and R13) arid the 


cfustering measurements from BOSS (White et al.|2012) 
to inform the properties of individuaf quasars for our 
simulation input. By using these two constraints, as well 
as a formalism for populating dark matter halos with 
quasars that we will outline below, we are able to select 
simulated quasar hosts that agree well with the latest ob¬ 
servational constraints. Specifically, matching the QLF 
means that we have an observationally accurate number 
of ionization sources, and matching the clustering mea¬ 
surements means our topology of reionization {e.g., the 
size and overlap of reionized regions) will be similar to 
the actual reionization process. The clustering can also 
have an effect on the spatial correlations present in the 
radiation field, which can affect the baryon acoustic os¬ 
cillation (BAG) measurement from the Lyman-a forest. 

This first paper of the series discusses the way in which 
we create sources for our simulations of helium reion¬ 
ization. In Sec. we describe our simulation strat¬ 
egy, and how we construct a quasar catalog from an 
Wbody halo catalog. In Sec. ^ we explain how we 
modify our quasar properties in order to match recent 
observations. In Sec. we explore implications of our 
findings for quasar populations. In Sec. we discuss 
implications for helium reionization. Finally, in Sec. 
we summarize our presentation and lay out future direc¬ 
tions. Throughout this work, we assume a ACDM cos¬ 
mology with rim = 0.27, Da = 0.73, Df, = 0.045, h = 0.7, 
as = 0.8, and Yne = 0.24. The se values are consist ent 
with the WMAP-Q year results ( Hinshaw et al.p013 ). 


2. MODELING QUASARS AS RADIATION 
SOURGES 

2.1. Radiation-hydrodynamic simulations 

When modeling helium reionization, we employ the 
RadHydro code, which includes Af-body, hydrodynamics, 
and radiative transfer calculations simultaneously. The 
code includes a particle mesh (PM) solver for gravity cal¬ 
culations, a fixed-grid Eulerian code for solving hydrody- 
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namics, and radiative transfer solved by performing ray¬ 
tracing. For more details on the hydrodynam ics portion 
of the simulation code, seejTrac & Pen (2004). For more 
details regarding the RadHy dro code and its app lic ation 


to hydrogen reionization, see Trac & Cen (2007) or Trac 


et al. (2008). 

i'he simulation strategy we employ for our exploration 
of helium reionization consists of two steps. First, a high- 
resolution A'^-body simulation is run for a given set of 
initial conditions. Halos are found on-the-fly using the 
friend-of-friends algorithm, and a corresponding catalog 
of spherical overde nsity halos are sav ed at even steps in 
cosmological time (Trac et al. 2015). Then, using the 
same initial conditions, a medium-resolution simulation 
using the RadHydro code is run. In order to provide 
accurate sources of ionizing photons for the radiative 
transfer calculations, it is necessary to convert the halo 
catalogs produced from the first simulation into quasar 
catalogs for the second simulation. Since the resolution 
of the RadHydro simulations is comparatively low (typi¬ 
cally a hydro grid unit is 10-100 /i“^kpc), the simulations 
are not able to accurately capture the subgrid, galaxy- 
level physics to include quasars directly. Thus, either 
a halo-level scaling relation or observational constraint 
is needed in order to create a physically reliable sam¬ 
ple. Rather than having to rely on scaling relations that 
require several steps to convert between halo mass and 
quasar luminosity, we use abundance matching to cal¬ 
culate luminosity as a function of mass, and then use 
observations to create a population with the proper char¬ 
acteristics. 

In order to calibrate the proper quasar properties 
to use, a suite of 10 IV-body P^M simulations with 
L = \ hr^ Gpc and 2048^ dark matter particles were 
run, which corresponds to a particle mass of rup = 8.72 x 
10® h~^MQ. The total volume is thus 10 {h~^ Gpc)^; the 
BOSS measurement of th e two-point correlation func¬ 
tion in White et al. (2012) has an effective volume of 9.8 
{h~^ Gpc)'^, so the volumes are comparable. Then halo 
finding was performed which produced the associated 
halo catalog snapshot every 20 Myr between 2 < z < 10. 
Since only comparatively massive halos serve as hosts 
for the bright quasars of interest, the simulations have 
a sufficient resolution to capture the required number of 
halos. 


2.2. Quasar light curves 

The first step in our model construction is to define 
the properties of individual quasars. The two most im¬ 
portant of these are the light curve (ie., L{t)) and the 
quasar lifetime. The most common model found in the 
literature for the light curve of quasars is the so-called 
lightbulb model, in which a quasar emits radiation at a 
constant luminosity for a lifetime tq before turning off. 
Though largely unphysical, this model has the conve¬ 
nience of being simple to implement in calculations. A 
further simplification is typically made in which it is as¬ 
sumed that tq is independent of luminosity, so that this 
quantity becomes a universal property. 

A more realistic model of the light curve is to assume 
an exponential form. This type of model can be moti¬ 
vated physically by noting that it corresponds to Edding¬ 
ton accretion onto the central SMBH. Several variations 


on this version include an exponential ramp-up to some 
peak luminosity followed by abrupt turn-off, a symmet¬ 
ric exponential about some peak luminosity, or an ex- 
pon ential ramp-up with a power-law fall-off in luminos - 
ity (jHopkins &: Hernquist 2009 McQuinn et al. |2009 ). 
While these models are more physically motivated, they 


are slightly more complicated. The approach we outline 
below is able to reproduce a given luminosity function 
for quasar light curves of this form. 

Specifically, we consider here two classes of quasar light 
curves: the “lightbulb” model and “exponential” model, 
defined as: 


Llh{t) — Tpeak0(t + ^g/2 — to)Q{tq/2 — t + tg), (1) 
Texp(t) — Apeak 6Xp( 1^0 (2) 

where 0(t) is the Heaviside theta function and tg is the 
time when the quasar reaches its peak luminosity Apeak- 
In the exponential case, the parameter r can be treated 
as a free parameter in a manner analogous to tq in the 
lightbulb case. Nevertheless, we re late t to tq, which we 
will describe in more detail in Sec. 12.41 
Another consideration is the quasar lifetime itself, 
which in general need not be a universal property of 
all quasars. We have parameterized quasar lifetime as 
a function of luminosity using a power-law form: 

where we vary the values of tg and 7 . We explore models 
in which 10^ < to < 10® yr, and -0.25 < 7 < 0.10. Pos¬ 
itive values of 7 imply that brighter quasars have longer 
lifetimes compared to dimmer ones, and 7 = 0 is the case 
of a universal lifetime for all quasars. 


2.3. Triggering rate 

We have discussed considerations for the individual 
quasars (ie., light curves and lifetimes), and we wish 
to connect them to the universal quasar population (ie., 
the QLF). In order to do so, we use the concept of a 
triggering rate h(Apeak 5 z), which dictates the differential 
number density of quasars that reach their peak luminos¬ 
ity Apeak as a function of luminosity and redshift per unit 
logarithmic luminosit y. Using the formalism outlined in 
Hopkins et al. (2006), we distinguish between the peak 
luminosity of a quasar Apeak and the instantaneous lumi¬ 
nosity at which it is measured for the construction of the 
QLF A, and relate the two with the triggering rate n. Es¬ 
sentially, the triggering rate must be convolved with the 
light curve of the quasars, since the measured luminos¬ 
ity function reflects a given quasar’s current luminosity 
A rather than its intrinsic peak luminosity Apeak- The 
result of this convolution is the observed QLF from the 
intrinsic triggering rate: 

(j){L,z)= J A (Apeak, z) d log Apeak - (4) 


As explained in Hopkins et al. (|2006 ), (/)(A) is the 
QLF (ie., the comoving number density of quasars 
per logarithmic bin in luminosity), and the quantity 
dt{L, Apeak)/dlog A is the amount of time that a quasar 
spends in a logarithmic luminosity bin. Essentially, the 
triggering rate can be thought of as analogous to the halo 
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mass function, though with the light curve convolution 
to account for changes in quasar brightness. In simple 
cases of the light curve the triggering rate can be solved 
for analytically: in the case of a lightbulb light curve, 
dt(L, Lpeak)/dlogL is a delta function at L = Lpeak) 
and so the triggering rate is proportional to the QLF: 

hlightbulb(I^i (5) 

tq 


In the case of an exponential light curve as defined in 
Eqn. ([^, we have 


Rexp(E, z') 


J_ d(/)(L,2:) 

2t d log L 


^ —-^peak 


( 6 ) 


activation of the quasar. Thus, we need to introduce a 
factor to account for the fact that not all halos host an 
active quasar. If we assume that the fraction of halos 
hosting an active quasar is universal (ie., independent 
of halo mass or quasar luminosity), we can express abun¬ 
dance matching for quasars, assuming a lightbulb light 
curve, as: 

(j)(> L) = /onnhalo(> M). (8) 

Expressed this way, /on is simply the fraction of halos of 
a mass M that host an active quasar. Alternatively, we 
could define this fraction in terms of the quasar lifetime: 

/o„(L,z) = ^, (9) 

tH{z) 


where the factor of 2 arises because a quasar will be ob¬ 
served at a luminosity L while its luminosity is increasing 
and then decreasing. In practice, the QLF is typically re¬ 
ported in magnitude units rather than luminosity. One 
common convention is to report the quasar’s absolute i- 
band magnitude at z = 2, Mi{z = 2). This quantity 
is then converted to the specific luminosity at 2500 A, 
units (erg s“^ Hz“^) by using Eqn. (4) of 


Richards et al. ( 2006| ): 


log 


10 


2500 A 
47rd^ 


- QA[M,{z = 2)+ 48.60 -f 2.5 logio(I + 2)], (7) 


where d = 10 pc = 3.08 x 10^® cm. To find the approx¬ 
imate bolometric luminosity, the relation of |Shen et al.| 
(2009) can be used to convert Mi{z = 2) to luminosity 
in erg s“^: 

Mi{z = 2) = 90 - 2.5 logio(A)- 


One should note that this relation is approximate, and 
depends on the assumed spectral energy distribution 
(SED) of the quasar. Eqn. (|^ is soluble for a few classes 
of light curves, such as the ones explored here. 


2.4. Abundance matching 

The technique of abundance matching has already been 


applied to populations of galaxies with { 

;reat success 
and has also 
0 ., Martini & 

{e.g., Simha et al.|2012 Hearin et al.|2013). 

been discussed in the context of quasars (e. 

Weinberg|2001 Porciani et al.|2004 Croton 

(2U0!)|). How- 


ever, we wish to extend the techniques mentioned above 
to include different quasar light curves and lifetimes. The 
methods we outline below are also fairly general, and can 
be extended to include semi-analytic models as well. We 
start with the Ansatz for abundance matching of galax¬ 
ies, namely that the most luminous galaxies are found 
in the most massive halos. This makes intuitive sense: 
more massive halos have more dark matter and baryonic 
matter to eventually convert to stars. Specifically, halo 
mass is highly correlated to the luminosity in the red 
bands, which shows the percentage of older stellar mass. 

For quasars, we have a similar situation where the most 
luminous quasars are found in the most massive halos. 
However, in this case the situation is slightly more com¬ 
plicated because quasars have a lifetime which is much 
shorter than the period from the halo’s formation to the 


lifetime (Martini & Weinberg 

2001 

), or the Hubble time 

(Conroy & White 

2l)I3|). we follow 

Conroy & White 

(2U13I and use the 

the exact choice fo 

Hubble time. As we shall see, though, 
r tniz) does not strongly affect the 


results. For the redshifts of interest, for a uniform value 
of tq = 30 Myr, this implies that /on ~ 0.1 — 1%. 

We can generalize the procedure of abundance match¬ 
ing to different light curves by using the triggering rate. 
In integral form, we can write abundance matching as 
equating the cumulative number of quasars above a par¬ 
ticular peak luminosity given by the triggering rate with 
the cumulative number of halos given by the halo mass 
function. The total number of halos which should host 
quasars within a time interval At is: 


[ [ fi{L*)d\ogL*dt 

JAt JL 

f /■~dnhaio(A*)dlogM*dP 


J At Jl d log M* d log L* dt 

M dnhaio(M* ) 

tn Jm dlogM* 


dlogL* dt 


( 10 ) 


■dlogM*. 


This form of our abundance matching equation becomes 
the central mechanism by which we are able to equate 
quasar luminosity with host halo mass. In this con¬ 
struction, we have implicitly used the mass-to-light ratio 
d log M /d log L to convert halo mass to quasar luminos¬ 
ity. Additionally, we have introduced the factor dP/dt 
to represent the probability that an individual halo will 
host a quasar. We have set this quantity to be equal 
to l/tn- Thus, for the case of a lightbulb light curve 
and a universal quasar lifetime, this formalism reduces 
to Eqn. (1^. Formally, this expression is an expansion 
oi hiL* about z th at is first-order accurate to At/tn 
(|Hopkins et al.||2006|). Thus, so long as the time-steps 
between deterniming the triggering rate are small com¬ 
pared to tn (defined either as the Hubble time or the 
halo lifetime, both several orders of magnitude longer 
than the typical quasar lifetime), this expression should 
reproduce the target QLF. 

In the exponential case, we are free to choose the pa¬ 
rameter T in any way that we like, as long as it is constant 
with respect to L (though it may vary with Apeak)- We 
have chosen r such that h(Apeak) is the same between 
the lightbulb and exponential cases for all luminosities. 
We accomplish this by equating Eqn. ([^ and Eqn. (|^, 
and solving for r in terms of tq. The expression involves 
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the ratio of the QLF and its derivative. This means that 
when we perform abundance matching, the same implicit 
mass-to-light ratio is used in the two cases. Since the 
halo mass function is the same between the two cases 
(due to the same population of halos being used), and 
the functional form of h is the same, we must have the 
same form of d log M /d log L. This has the advantage of 
allowing us to apply certain intuition from the lightbulb 
case to the less straightforward exponential case. The 
downside to this approach is that when exploring the 
parameter space of quasar lifetimes tg in the lightbulb 
case, it is not immediately obvious how this translates 
to the exponential time constant r, since we effectively 
have different values of tg for different luminosities. For 
instance, even in cases where tg is independent of lumi¬ 
nosity, r still changes as a function of L. However, the 
benefits of being able to interpret the results of the expo¬ 
nential case using the intuition provided by the lightbulb 
case outweigh the downsides of not exploring parameters 
in T directly. 

The general procedure is as follows. 


1. The halo mass found from the halo catalog at red- 
shift Zcat is read and converted to an expected 
number density in a particular cosmol ogy using the 
univer sal mass function described in ITinker et al.l 
(20081. The fitted form of the mass function is 
used rather than the empirical one from the cat¬ 
alog in order to decrease the variation in number 
density at the high-mass end, since these quasars 
are disproportionately important for the reioniza¬ 
tion process. 


2. Using Eqn. (10), the halo number density is con¬ 


verted to an expected quasar number density using 
a specified QLF. 


3. The quasar magnitudes are binned into equal inter¬ 
vals in magnitude AM, such that the expected trig¬ 
gering rate h(M, Zcat) ■ • ■ h(M-|-AM, Zcat) is found, 
which is converted from a number density to a total 
number N (M) using the volume of the simulations. 


and one of them is hosting an active quasar, the de- 
scendent halo inherits the active quasar. If a single 
active progenitor halo splits to form two descendent 
halos, the larger halo retains the quasar. In the case 
of a merger between two active quasar halos, only 
the larger quasar survives. These cases represent a 
comparatively few number of instances of our total 
population evolution, and do not strongly influence 
our conclusions. 

2.5. The quasar luminosity function 

Throughout this work, we use a series of QLFs as de¬ 
termined at different epochs. For relatively low-redshift 
(2 ^ ^ 3), we use the QLFs as determined by R13 from 

the BOSS survey, specifically the high-z stripe 82 sample 
(S82) form which includes luminosity evolution and den¬ 
sity evolution (LEDE). Above a redshift of 3, the QLF 
has been measured at z ^ 3.2 and z ^ 4 by M12 us¬ 
ing data from COSMOS. At z ^ 5, the QLF has been 
measured by M13 using data from the SDSS.j^ Although 
these works use slightly different values for cosmological 
parameters from the ones assumed here, the impact on 
the reported quantities is minimal. 

In order to span the different epochs over which the 
luminosity function has been measured, it is necessary 
to combine the different data sets. All of the data sets 
fit to a double power law form of the QLF, written as: 

(j)^ 

— j^qo.4(i+q)(M-M*) _|_ 200-4(1+/3)(M-M*) ’ 

where $ is the comoving number density of quasars of 
magnitude M per unit magnitude, 4>* is the normaliza¬ 
tion of the QLF, a is the faint-end slope of the luminos¬ 
ity function, /3 is the steep-end slope (which is reversed 
from the parameterizations of M12), and M* is the so- 
called break magnitude where the luminosity function 
transitions from the faint-end to the steep-end. In most 
formulations at high-redshift, redshift evolution is incor¬ 
porated by a change in (ft*, M*, or both, that is linear in 
redshift. For the data from R13, the evolution is given 
by the equations: 


4. Within each magnitude bin, each quasar is assumed 
to have an equal probability of becoming active. 
Each quasar candidate is randomly turned on with 
probability l/Nun{M). 

5. To ensure that the volume self-consistently fol¬ 
lows the merging of the underlying host halos, the 
quasars are propagated forward using a halo merger 
tree. By design, the halo catalog snapshots are 
made at times that are shorter than the expected 
lifetimes of the quasars. This approach allows for 
halos hosting quasars to be tracked throughout the 
simulation. In most cases, an active quasar from 
time step * — 1 in a progenitor halo passes to the 
single descendent halo at time step i. Additionally, 
this halo hosting an active quasar is not eligible 
to host a new quasar. This approach covers the 
majority of halos for the majority of time steps. 
However, there are several special cases related to 
merger events worth discussing. Specihcally, when 
two progenitor halos merge into a single descendent 


logio = logio ^0 +ki{z- 2.2), (12) 

M*{z) = M* + k2{z-2.2). (13) 

For the data in M13, there is linear evolution in logj^g f* 
as well, given as: 

logio(/)*(z) =logio«!>5 + fci(z- 6 ). (14) 

To combine the R13, M12, and M13 data sets into a 
single set of quantities, we first assume that the results 
from R13 are accurate for redshifts z < 3.5. This is the 
nominal limit of the LEDE fits, and thongh there are 
small differences between the fit QLE and the binned 

^ Additional l y, the QLF at 2 : ^ 4 has also been measured by 
|Glikman et al.| (j2011ll andllked a et al1l|201 1 ||. As noted in M12, 
the normalization ot the QLh otllkeda et al.lJzOlT ^ is comparable, 
whereas the normalization of|Glikman et aT| (|2Uil|| is larger than 
the others by a factor of ~ 4. NLL'Z notes that the difference can 
be caused by contamination of the faintest-magnitude bins from 
dwarf stars and high-redshift galaxies. In the following analysis, 
we use the results from M12. 

^ A n upper limit for the QLF at 2 : ~ 5 was found by|lkeda et al.| 
< |2012[ |, which is consistent with the results of M13. 
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Table 1 

A list of the QLF parameters of the Datasets Incorporated. 


Dataset 

z 


M*b 

ki^ 

k2 

a. 

/3 

R13 

M12 

M12‘i 

M13® 

M13 

M13 

2.2-3.5 

3.2 

4 

5 

5 

5 

-5.93l°;°? 

-6.58l°;?® 

—7 

—8 47+'^-20 
-0.24 
•j ^q-|-0.30 
< •Oo_o 25 
■7 Qq-|-0.03 
' •^'^-0.03 

—26 

zv.oi _o.o2 
-27.03 ±0.68 
-27.13 ± 2.99 
-28.701°;"^ 
-27.34l°;t° 
-27.88 

-0.689l°;°27 

-0.47 

-0.47 

-0.47 

-0.809lO;“3 

-*-•^^-0.03 
-1.73 ±0.11 

-1.72 ±0.28 

— 2 

^•'^'^-0.14 

-1.50 

-1.80 

0 cri+0.09 

-2.98 ±0.21 

-2.6 ±0.63 

-4.00 

q -1 Q-|-0.28 

'^•-^^-0.41 

-3.26 


^ (p* has units of Mpc ^ mag ^. 
b M* = Mi{z = 2) = Mi450 - 1.486. 

^ k\ and /c 2 are defined for models with redshift evolution in Eqns. l |12|14[ |. 

^ The authors of M12 provide a value for 0 q at 2 ~ 4 where the reported error is greater than the value itself. Since 
this value must be positive, the resulting lower-bound is unphysical. We reproduce the value and upper-bound here 
for completeness, but do not include this value directly when determining the values of the QLF. See Appendix [A| 
for further details. 

® In M13, the authors provide three fits, each with at least one parameter held constant. Values without error 
ranges indicated correspond to the parameters held fixed for a particular fit. 


data, overall the fits are excellent. To incorporate the 
results at higher redshifts, we cast the four parameters 
of the QLF M*, a, and /3) as quantities that have 
linear evolution in redshift. We define these parameters 
as: 


logio = logio<i!>S + ci( 2 :- 3.5), 
M*(z) =M* + C2 (z-3.5), 
a{z) = ao + C 3 (z - 3.5), 

/3(z) = /3o + C4(z - 3.5). 


Table 2 

A List of the Parameters Used in Eqns 
Based on the Data Listed in Table 


. |^5a)15d] l 



Parameter 

Fiducial Value 

Parameter Range 

(15a) 

logic 

Cl 

-6.82 

-0.790 

[-1.10,-0.536] 

(15b) 

M* 

C2 

-27.6 

-0.238 

[-0.716,0.170] 

(15c) 

QO 

-1.29 

C3 

-0.324 

[-0.493,-0.140] 

(15d) 

^0 

-3.51 


C4 

0.0333 

[-0.327,0.260] 


These parameterizations are applied to redshifts where 
z > 3.5. The constant values are defined to be equal 
to the values of R13 at z = 3.5, and the values for the 
slopes (C 1 -C 4 ) are allowed to take on a range of values. 
The range is generally chosen such that the values for the 
different parameters brackets the range of best-fit values 
provided by the highest redshift (M13) data. The fidu¬ 
cial values for the slopes are taken to be ones that reason¬ 
ably reproduce the high-redshift measurements. Table 
shows the fiducial values for the slopes, as well as the 
range of values for the parameters at z ~ 5 used in the 
parameter space exploration in Sec. For a complete 
discussion on selecting the parameters for the QLF, see 
Appendix [A] 

Table [^lists the parameters that we include from the 
measurements of R13, M12, and M13. The values from 
M12 are not included in the fitting procedure directly, 
and serve primarily as a consistency check due to their 
comparatively large error bars. The parameters from 
Ml3 are determined at z ~ 5, and the ones from M12 
are determined at z ^ 4 and z ^ 3.2. Note that the 
authors of M13 provide three independent fits to their 
data, which are all incorporated into the final QLF pa¬ 
rameterization. (See Appendix [A] for more details.) For 
the measurements from R13, whose fiducial LEDE model 
includes redshift evolution in </)* and M*, the model is 
valid over a range of redshift, from 2.2 < z < 3.5. Eor 
the purposes of generating our quasar catalogs, we are 
interested in exploring the QLF until z = 2. For the 
sake of simplicity, we simply extend the LEDE model 
from R13 to this redshift. Although the LEDE fit is 


Note. — These provide a fit to the lumi¬ 
nosity function through redshift, and ensure that 
the abundance of quasars matches observations 
as nearly as possible. For additional details on 
the parameters and the fitting procedure, see Ap¬ 
pendix^ 


ostensibly not valid below z = 2 . 2 , we expect helium 
reionization to be largely finished by this redshift, and 
so the precise form of the QLF at z ~ 2 is not of fun¬ 
damental importance to our study. Also, for the value 
of M*, it is necessary to conv ert t o a single magnitude 
system. As explained in Sec. 2.3 we use Mi{z = 2 ), 
the absolute t-band magnitude at z = 2. The QLFs of 
M12 and M13 use M 1450 , which i s related to Mjjz = 2) 


by Mjjz = 2) = Mi 45 o — 1.486 (Richards et al. 2006 


Ross et al. 2013 Appendix B). Note that this conver- 


sion assumes that the quasar SED follows a power-law 
with an effective spectral index of a = 0.5 (using the 
convention that fviy) oc Modifying the spectral 

index a changes the magnitude conversion, so care must 
be taken when converting between magnitude systems. 
See Appendix [A| for further discussion. 

Eigure shows the combined QLF from R13, M12, 
and M13 ^hich at this epoch is essentially that of R13), 
as well as two different quasar models at z ~ 2.4. We 
can see that there is generally very good agreement be¬ 
tween the constructed quasar catalog and the target lu¬ 
minosity function, as should be expected. The differences 
between the constructed catalogs and target luminosity 
function are typically on average 5%, which is compa- 
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Figure 1. A comparison of the composite quasar l uminosity func¬ 
tions from the SDSS+COSMOS measuremen ts l|Masters et al.| 
|2012[ [McGreer et al.]|2013| |Ross et al.||201^ to our abundance 
matcliing mettiod, plotted with Foisson error bars. The two differ¬ 
ent quasar models (defined in Table are offset from each other 
for visual clarity. The agreement is excellent for comparatively dim 
quasars which are more common, but there is some discrepancy for 
bright objects. The reason for this disagreement is primarily due 
to Poisson noise, since these objects are rare even for the large (1 
(h~^ Gpc)^) simulation volume. At low luminosity in the expo¬ 
nential case, the completion limits of d ark matter halo hosts at this 
mass become noticeable. See Sec. El for further discussion. 

rable to or smaller than the uncertainties in the luminos- 
ity function itself at these redshifts. At high luminosities 
{Mi < —28), though, there are some comparatively large 
differences that can arise between the predicted and em¬ 
pirical luminosity functions. This deviation is largely 
due to Poisson shot-noise introduced by the rarity of 
the objects. For objects in this luminosity range, there 
are typically only a few objects ( 0 ( 10 )) in the entire 1 
(h~^ Gpc)^ volume. At the dim end of the QLF, there 
can be insufficient halos of a particular mass given the 
mass resolution of our simulation. The minimum halo 
mass is Mhaio.min = 4.36 x 10^^ M^Mq. Since quasars 
with Mi < — 25 are most important for this study, this 
does not affect our results significantly. 

Throughout most of the following analysis, we focus 
our attention on several models in particular, parame¬ 
terized in terms of tg and 7 as in Eqn. ® . The first four 
of these models have particularly goooagreement with 
the BOSS measurements. The last two are included to 
demonstrate how the clustering signal changes as a func¬ 
tion of to for a fixed value of 7 : models LI, L3, and L4 
all have the same 7 value. We summarize these models 
in Table m 


3. CLUSTERING MEASUREMENTS 
3.1. Two-point Correlation Function 


By construction, our method matches the input QLF 
at all redshifts, regardless of the individual properties of 
the underlying quasar population. However, we are not 
guaranteed to match the observed clustering of quasars. 
Changing the implicit mass-to-light ratio of Eqn. (10 1 
through changing the quasar lifetimes will affect how na- 
los are populated with quasars. In general, longer quasar 
lifetimes lead to quasars of the same luminosity being 
matched into hosts of larger masses. Since their hosts 
are more biased, this leads to quasars of the same lu- 


Table 3 

A List of the Parameters of Some Quasar Models 
Considered. 


Model Name 

Light Curve 

logic (A/yr)*" 

7 

LI 

Lightbulb 

7.75 

0 

L2 

Lightbulb 

8.25 

-0.125 

El 

Exponential 

7.25 

0 

E2 

Exponential 

7.75 

-0.15 

L3 

Lightbulb 

7 

0 

L4 

Lightbulb 

8.5 

0 


^ £0 ^-iid 7 as defined in Eqn. 


minosity showing a larger clustering signal. This is true 
at all luminosities. We want to match the clustering be¬ 
cause it can affect the topology of reionization. There 
can also be spatial correlations present in the radiation 
field as a result of reionization, which are important for 
making me asurements of the BAO from the Ly man-g 
forest {e.g., White et al.||2010 Slosar et al.|2013 ). 

Here, we explore how to include clustering measure¬ 
ments from the two-point correlation function in our 
quasar catalog. Recent results from the BOSS survey for 
the clustering of quasars in the redsh ift range of interest 
are presented in White et al. (2012 1 . The above work 
examines the clustering signal of quasars in both 2D- 
projected and 3D-redshift-space correlation functions at 
intermediate scales (3 < s < 25 /i“^Mpc). The authors 
also introduce luminosity cuts to make the results more 
robust. For the purposes of this comparison, we consider 
their selection for which they imposed luminosity cuts 
on both the bright and faint ends, so that only objects 
with —25 > Mi > —27 were considered across the entire 
redshift range (Sample 4 as defined by the authors). For 
a fair comparison, we impose similar cuts on our object 
selection. We also examine the redshift evolution of the 
results, and compare against the high-z/low- 2 ; samples 
(Samples 5 and 6 ) as well. See Appendix for further 
discussion of these different redshift samples. 

We explore the parameter space of available quasar 
models by examining the lightbulb and exponential light 
curves defined in Eqns. ([^ and ([^, as well as luminosity- 
dependent quasar lifetimes dehned in Eqn. ([^, param¬ 
eterized by tg and 7 . Eor each combination of param¬ 
eters, we construct a quasar catalog in the manner de¬ 
scribed abovej^ Then, we extract from this catalog all 
objects that satisfy the magnitude constraints at the cen¬ 
tral redshift of the survey z = 2.39. This redshift repre¬ 
sents the average redshift of quasars chosen in the BOSS 
sample; the actual quasar objects s pan in redshift from 
2.2 < z < 2.8. However, as noted in White et al. (2012), 
the redshift evolution of the signal is weak. Thus, ex- 


® There are several extreme models where the number of objects 
is significantly fewer than the number predicted by the quasar lu¬ 
minosity function. This is not a failure of our methodology, but 
rather instances of there being too few halo objects of a given mass 
to host quasar objects. In essence, /on is so small that we reach the 
resolution limits of the simulation. In these cases, we add particles 
from a second-order Lagrangian perturbation theory (2LPT) sim¬ 
ulation of the same initial conditions at the same redshift in order 
to define a set of “random” particles that are still representative 
of the underlying matter distribution. We randomly sample from 
these particles in order to fill out the catalog to the expected num¬ 
ber. This ensures that we do not measure a statistically significant 
clustering measurement when the catalog is clearly unphysical. 
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Figure 2. The correlation matrix for the LI model. Note how 
the matrix is dominated by the diagonal entries, which is to be 
expected for shot-noise dominated measurements. The small off- 
diagonal terms suggest that the covariance matrix has converged 
numerically, and should be stable when inverting. This type of 
structure is seen in all models considered. 


tracting objects from our quasar catalogs at a single red- 
shift rather than a range should have little effect on our 
overall conclusions. We measure the monopole of the 
two-point correlation function using the “natural esti¬ 
mator” 


as) 


{DDjs)) 

{RR{s)) 


(16) 


where {DD{s)) is the average number of quasar pairs 
from the quasar catalog separated by a real-space dis¬ 
tance of [s —As/2, s-l-As/2], and {RR{s)) is the number 
of pairs of points at the same separation drawn from a 
distribution with Poisson noise. 


3.2. Calculating values 

In order to quantify the statistical uncertainty in our 
catalog, we ran a suite of 10 A-body simulations with 
different initial conditions. We then performed our abun¬ 
dance matching procedure on each of the different sim¬ 
ulations, including several realizations for each volume. 
Since our abundance matching procedure stochastically 
determines which halos should be hosting active quasars 
at a given time step, we create several quasar catalogs 
for each individual halo catalog, using a different initial 
random seed (three realizations per volume for these re¬ 
sults). Additionally, we have augmented the effective 
number of samples by including redshift space distor¬ 
tions along the different principal axes of the simulation. 
This strategy gives us a total of 90 samples for which to 
measure the clustering signal. The best estimate for the 
correlation function ^(s) for a given radial bin Si is given 
by averaging over all of the individual estimates ^k- 

1 ^ 

a^^) = JJT.a{s^) ( 17 ) 

' fc=i 


We then estimate the covariance between the radial bins 
by computing the entries of the covariance matrix Cg. 


We compute the entries of the covariance matrix as (Ze- 



Figu re 3. The quasar two-point correlation function from |White| 
|et al.| l |Ml2|| , collared to several models whose parameters are de¬ 
scribed m 'iable[3| All measurements were made at the same values 
of s, but are offset from each other for visual clarity. The shaded 
error regions on the measurements from BOSS are the reported Icr 
error bars, and the error bars on the models are the square root 
of the diagonal elements of the covariance matrix. Note that for 
the same value of 7, increasing £0 leads to a larger clustering signal 
(compare L3, LI, and L4 in order of increasing £0). See the text 
for additional details. 


havi et ar]|2QQ5 1: 


N 


= J^'^{a{si)-as^))ia{sj)-asj))■ (is) 


k=l 


The correlation matrix entries for our model LI is plot¬ 
ted in Figure Notice that the diagonal entries dom¬ 
inate, which means that the bins are mostly indepen¬ 
dent of each other and dominate d by shot-noise (Valageas 
et al.||2011 White et al. 20121. Implicitly, the samples 
have been treated as being independent, and this is al¬ 
most surely not the case. Although the 10 volumes as 
a whole can be treated as being statistically indepen¬ 
dent, the different realizations based on the same halo 
catalog are likely correlated. Further, the projections 
of peculiar velocities along different axes for the same 
realization are also likely to produce correlated results. 
However, producing a sufficient number of independent 
realizations to decrease the noise in the covariance ma¬ 
trix is computationally infeasible. Further, the variance 
in the clustering signal among quasar catalog realizations 
for a given (to)7) pah is comparable to small displace¬ 
ments in the to“7 parameter space, so it is necessary to 
include this source of uncertainty. Since we are interested 
only in finding models that are consistent with the BOSS 
measurements which have their own set of observational 
uncertainties, we feel that this approach produces suffi¬ 
ciently accurate results. 

Once the entries of the covariance matrix have been 
computed, the difference vector S{si) = Cmodei(si) — 
^BOSs('Si) is calculated. The correlation function ^boss 
is fit to a power law: ■Cboss(s) = (s/sq)^, where the 
authors have fixed the value of /3 = — 2. In order to in¬ 
vestigate the impact this choice has on the conclusions, 
we performed fits on the correlation function measured 
from our quasar catalogs using two different parameter- 
izations: one where the best-fit value of sq was found 
when fixing /3 = —2, and another where the value of 
So and /3 were both fit. In the length scales used for 
our analysis (3 < s < 25), the deviation of f3 from the 
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Table 4 

A List of the Best-fit Parameters for Our Quasar Model as a 
Function of Redshift. 


Redshift Selection 

■^eff 

Light Curve 

T * a.b 
''"efl 

f* C 

High-2d 

2.51 

Lightbulb 

Exp 

12.92 

12.40 

7.62 

7.14 

Fiducial 

2.39 

Lightbulb 

Exp 

13.29 

13.05 

7.77 

7.18 

Low- 2 : 

2.28 

Lightbulb 

Exp 

13.17 

13.15 

7.84 

7.29 


^ Leff and tgff as defined in Eqns. ( |20|2I| |. 

^ ^eff “ ^'^Sloi^eS/L q) 
ieff = logl0(ieff/yi-) 

The high -2 and low -2 samples examine the evolution of these pa¬ 
rameters with redshift. See Appendix^for further discussion. 


fiducial value of —2 was small, typically less than 5%. 
Furthermore, the values for sq were also largely similar 
between a fixed slope or a varying one, with deviations 
typically less than 1%. Thus, the choice to set (3 = —2 
does not strongly bias the results presented here, or the 
values reported in ^boss- 

When comparing one of the quasar models with the 
BOSS results, the value of the model is then given 
by: 

(19) 


To define the model that fits the BOSS observations best, 
we want to minimize the value of the model. A two- 
dimensional space in to and 7 is constructed for both of 
the light curves, and this space is explore d using regu¬ 
lar gr id points. Following the analysis of [White et al.| 
( 2012 ), a distribution with nine degrees of freedom is 
assumed. Using this distribution, the x^ value for a par¬ 
ticular model is converted to a confidence interval. An 
equivalent na value is computed based on the confidence 
interval (Icr if the enclosed probability is 0.683, 2cr if it 
is 0.955, etc.). This statistic demonstrates how “consis¬ 
tent” a particular model is with the BOSS observations. 

Figure shows the clustering measurements for sev¬ 
eral of our well-htting models compared to the BOSS 
measurements. The values of these models are given in 
Table In general, as tp increases at a hxed value of 7 , 
the clustering signal increases as well. Compare specif¬ 
ically the L3, LI, and L4 models, which have the same 
value of 7 but have respectively increasing values of to- 
Mathem atic ally, this behavior can be seen from the form 
of Eqn. (10): for the same luminosity and mass functions 
but a larger value of /on oc tq, quasars of the same lu¬ 
minosity will shift to more massive host halos. Since the 
clustering signal increases with the mass, it follows that 
increasing to will increase the clustering signal. For sim¬ 
ilar reasons, increasing values of 7 for constant values of 
to are also associated with a stronger clustering signal, 
since this also effectively increases the quasar lifetime tq. 


3.3. Characteristic luminosity and lifetime 

Figure shows the values in the two-dimensional 
parameter space to and 7 , as defined by Eqn. ([^, for 
the different light curves. The region of good agreement 
between the BOSS measurements and our models takes 



Figure 4. A comparison of the parameter space exploration in 
terms of the parameters Iq and 7 from Egn. Both the parame¬ 
ter space for the lightbulb model (Eqn. (m) and exponential model 
(Eqn. are shown. The dashed lines represent the best linear 
fits to me data for a particular light curve. The class of models 
that are consistent with the BOSS measurements at Icr and 2a cor¬ 
respond to the darkly and lightly shaded regions. In general, we 
find that for the exponential model, shorter lifetimes are preferred 
(smaller values of to for tfio same 7). Since we abundance match 
against the quasar’s peak luminosity, and the quasar spends com¬ 
paratively little time at or near the peak luminosity, we effectively 
increase the clustering signal for lower luminosity quasars. 


on a linear relationship between log;^o(^o) and 7 . Such a 
relationship can be parameterized as: 


logio(^o/yr) = logio(teff/yr) + To7- (20) 

The parameters teff and Lq can be thought of as a char¬ 
acteristic timescale and a characteristic luminosity, re¬ 
spectively. From the functional form of our power-law 
for quasar lifetime in Eqn. ([^, Lq can be interpreted 
as changing the normalization luminosity. This is the 
luminosity at which all models have the same lifetime, 
regardless of the value of 7 . In other words, the charac¬ 
teristic luminosity of the power law becomes: 


^ogioiLefi/L q) — 10 — Lq, (21) 


where Lq is defined in Eqn. (20). The parameter tes is 
the characteristic time because all models have this same 
lifetime at the luminosity Teff- 
For the lightbulb model, the best-fit values are 
logio(teff/yr) = 7.76 and logio(Aeff/A 0 ) = 13.29. (See 
Table E for evolution of these parameters with redshift.) 
The characteristic luminosity inferred from this value is 
Lgff = which has a corresponding magnitude 

of Mi = —27.2. This value is not surprising, given that 
quasars were selected for the clustering measurements 
near this magnitude range. More interesting is the value 
of logio(<eff/yi') = 7.77, which gives a characteristic life¬ 


time of = 59 Myr. This is a quasar lifetime that 


is slig 
eratui 

btly longer than those tyi 

pically quoted in the lit- 

e IYu & Tremaine [2002 

Porciani et al.| 

2004 Yu 

& Lu 



2004 Conroy & White 2013p, which are 

closer to 


the balpeter e-tolding time scale or shorter (~45 Myr for 
a quasar accreting at Eddington luminosity and a mass 
conversion efficiency of e = 0.1). Although teS is slightly 
higher than these values, it is within a factor of 2 . 

In the exponential model, the best-fit values for t^s 
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Figure 5. The range of quasar lifetimes for models with the best- 
fit values of to as a fu nction of 7, based on the linear relationship 
extracted in Eqn. ( |20| l. The solid regions show the span in quasar 
lifetimes, while the clashed line shows the median value within a 
given model. Note that the median value is fairly constant across 
all quasar lifetimes. Thus we are able to characterize a quasar 
model reasonably well using the characteristic lifetime. The com¬ 
parative large spread in quasar lifetime in the exponential case is 
due to our method of selecting fg, rather than reflecting a truly 
large spread in the data. See the text for additional details. 


and Leff defined in Eqn (20l are logiQ(ieff/yr) = 7.18 
and log;^Q(Leff/ 7 v 0 ) = 13.05. This luminosity implies a 
slightly dimmer char acte ristic luminosity (M, = —26.6). 
As discussed in Sec. |2.4[ there is not a single r for all 
quasars for a given value of tes'. L m L* quasars have 
r « teff, with brighter quasars having r > teff- However, 
the difference between r and t^g does not differ by more 
than a factor of 2 in either direction, and so to a good ap¬ 
proximation r ^ teffi especially for the luminosity range 
used to match the clustering measurements. Compared 
to the lightbulb case, the quasars with an exponential 
light curve have a shorter characteristic lifetime of 15.1 
Myr. The characteristic lifetime is smaller for the expo¬ 
nential than in the lightbulb case because quasars do not 
shut off entirely after a single lifetime, so the time that 
a quasar is “bright enough” to be included within the 
luminosity cuts is longer than its lifetime t^g. This life¬ 
time is about a third of the Salpeter e-folding time scale, 
which implies that if quasar light curves are roughly ex¬ 
ponential, the combination of the measured QLF and the 
clustering measurements favors quasars that either radi¬ 
ate at luminosities dimmer than their Eddington ratio 
(L/Ledd < 1 ), have a mass-conversion efficiency that is 
less that the fiducial value (e < 0.1), or both. Unfor¬ 
tunately, since our model does not track the underlying 
physics present, we are not able to distinguish between 
these two cases. 

The reason for the different best-fit values between the 
two models can be understood as follows. By construc¬ 
tion, we have hxed the lifetime of the exponential quasars 
such that their peak luminosity-to-mass ratio is the same 
as in the case of t he lightbulb for a given choice of to and 
7 . (See Sec. 2.4 for more details.) However, the mass- 
to-light ratios tor the two light curves are significantly 
different. This is due to the fact that the observed lu¬ 
minosity for an exponential quasar can be much smaller 
than its peak luminosity. A particular luminosity range is 


selected for the clustering measurements, but the cluster¬ 
ing of these quasars is tied to their peak luminosity rather 
than the observed one. Thus, quasars will tend to have 
higher clustering at a given luminosity in the exponential 
case compared to the lightbulb, since they spend compar¬ 
atively little time at or near their peak luminosity. This 
luminosity selection includes quasars with a higher peak 
luminosity than the chosen range (and thus a higher clus¬ 
tering signal), so we must also include quasars that have 
lower mass hosts to match the average clustering signal. 
This means that there is a larger spread in host mass 
compared to the lightbulb case. This behavior explains 
why the characteristic luminosity is slightly smaller for 
the exponential model compared to the lightbulb: there 
is an increased number of low-luminosity quasars occu¬ 
pying high-mass hosts. 

Figure [^shows the range of quasar lifetimes as a func¬ 
tion of m^el parameter 7 . The quasar lifetime is broadly 
similar across different model choices. The exponential 
model has a lower overall value due to the effect discussed 
above, i.e., that quasars from a higher peak luminosity 
will be included in the sample, bringing along a higher 
clustering signal. Since this is true for nearly all the 
quasars in the sample, there is an overall decrease in the 
selected lifetime of quasars. The large difference in the 
span of quasar lifetimes is due to the way that we have 
defined the quas ar li fetime in the exponential model. As 
discussed in Sec. |2.4[ the exponential lifetime r is selected 
such that the same relationship between host mass and 
quasar peak luminosity exists in the exponential case as 
in the lightbulb case. Even in a model where for the 
lightbulb tg is independent of L (i.e., when 7 = 0), the 
exponential model parameter t does have luminosity de¬ 
pendence. In general, quasars with luminosities above 
L* will have a lifetime longer than an equivalent lumi¬ 
nosity in the lightbulb case for the same choice of to and 
7 in Eqn. (1^, and those with low luminosities will have 
a shorter liiCTime. This choice for our model leads to the 
spread in lifetimes of a factor of ~5, as seen in the case of 
7 = 0. For models in which 7 > 0, there is a widening in 
the range of values. This is due to the fact that brighter 
quasars live longer than dimmer ones. Since our choice 
of quasar lifetime already enforces this difference, these 
models see an increased effect. Conversely, for 7 < 0, 
there are competing effects between brighter quasars be¬ 
ing less long-lived due to the choice of 7 , but still simul¬ 
taneously living longer than their lightbulb counterparts 
due to the choice of tg. The latter effect wins out, and 
these quasars end up having a signihcantly larger spread 
than in the lightbulb case. Note that the contours in this 
figure are smooth compared to Fig. because these are 
results lying along the best-fit line, and the figure shows 
the range in values rather than a single number (ie., the 
value) that fluctuates as a function of position in to 
and 7 . 


4. DISCUSSION 
4.1. Mass-to-light Ratio 

The combination of the QLF and clustering measure¬ 
ments produces an important set of constraints on the 
space of potential quasar models. Here we investigate the 
implications of these models. One important implication 
is the mass of a typical halo for a given quasar luminos- 
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Figure 6 . A comparison of the mass-to-light ratio between the 
different quasar models. We have computed this ratio for each 
model in our parameter space, and weighted their contribution by 
their corresponding value. The lines show the median value by 
weight and the shaded regions show ±1(T and 2cr. Note that for 
a given luminosity, a quasar in the exponential model is found in 
a halo with a smaller mass. This is due to the fact that quasars 
with a peak luminosity significantly greater than the observed one 
are included in the luminosity range selected for the clustering 
measurements. These hosts have a higher clustering signal than 
quasars with a peak luminosity in the luminosity selection. This 
means we must also select lower-mass objects as well. See the text 
for additional details. 


ity. It is trivial to predict this for the case of a lightbulb 
model, but less straightforward for the case of the expo¬ 
nential model. Here the peak luminosity Lpeak is used to 
define the mass-to-light ratio, since this is the quantity 
used in our abundance matching approach. This ratio 


pared with results of previous analysis {e.q.. 

Martini & 

WeinbergJ|2001 

Shen et al.||2007 

White et al. 

2012 ). 


of the mass of the host halo for the lightbulb and expo¬ 
nential models for all combinations of to and 7 - This plot 
shows the mass-to-light ratio of the entire catalog. The 
weight assigned to the luminosity as a function of mass 
LpeakiM) for a particular model i is given by a likeli¬ 
hood. We also find the ± 1 and 2 a values that enclose 
68 % and 95% of the likelihood. 

Figure shows the mass range as a function of the 
model parameter 7 . Note that the range is essentially 
constant with respect to 7 . By averaging the median 
mass across all values of 7 , a characteristic mass for the 
two models can be defined. This characteristic mass is 
2.5 X 10^^ for the lightbulb model, and 2.3 x 10^^ 

h~^MQ for the exponential model. These values are 
broadly consistent with pr evious studies of quasar clus- 
tering meas urements (e .g., Porciani et al. |2004 Groom 


et al. 1120051 ILidz et al.|[Ml)D[ |Porciani fc lNorberg||2Ul)6[ 


White et al. 2012 |). Since in all models the same clus¬ 
tering signal of the quasars is being selected, there is 
an implicit requirement for the hosts to lie within a cer¬ 
tain mass range. Additionally, the mass range in the 
exponential case is significantly larger than that of the 
lightbulb model. This is again related to the fact that 
due to the light curve, quasars with a higher clustering 
signal are included within the luminosity sample, and 
so there must also be lower-mass hosts included as well 



- Lightbulb 

- Exp 

^°-^0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 

7 

Figure 7. The selected mass range for our different models. For 
the lightbulb case, there is a relatively narrow range in mass. Since 
the clustering of quasars is fixed as a function of luminosity, there 
is a very tight relationship between observed luminosity and under¬ 
lying host mass. In the exponential case, the range is much more 
extended. Since the mass-to-light ratio is fixed to be the same for 
the peak luminosity, the evolution within the model means that 
there will be quasars with higher peak luminosity (and higher-mass 
hosts) selected by the evolution. 


to balance the average clustering strength. There is a 
significantly larger spread above the median mass than 
below. The reason for this asymmetry is due to the dif¬ 
ference in number density: since the high-mass objects 
are rarer, a comparatively smaller range in low-mass ha¬ 
los is necessary to make the c luste ring signal equivalent 
to the lightbulb case. See Sec. |3.3| for further discussion. 

In the case of the lightbulb model, the halo mass 
that corresponds to the selected luminosity range of 
quasars is relatively tightly constrained. For the mod¬ 
els that agree with the BOSS measurements at Icr, 
the average halo mass ranges from 1.35 x 10^^ to 
4.93 X 10^^ h~^MQ for hosts of quasars within the mag¬ 
nitude cutoff. For the exponential model, there is a much 
larger range in halo mass: for the collection of models 
that agree at Icr, the mass ranges between 6.69 x 10^^ 
and 5.85 x 10^^ hr^MQ, almost an entire order of magni¬ 
tude (compared to about half an order of magnitude for 
the lightbulb model). Also note that the mass range is 
much larger than in the lightbulb case. This fact can be 
explained by noting that there is evolution in the mass- 
to-luminosity ratio during the lifetime of the quasar. Fur¬ 
ther, the e-folding time for these models is comparatively 
long, with typical values being r « 40 Myr. This means 
that there are high-mass hosts included in the sample 
of quasars chosen for the clustering measurements whose 
quasars are not at their peak luminosity. Since these 
hosts have a bias larger than the value preferred by the 
BOSS measurements, this sample must necessarily in¬ 
clude hosts which have a smaller clustering signal, so 
that on average, the total bias agrees with BOSS. 

There is a systematic shift upward in the mass of the 
exponential case compared to the lightbulb. This shift is 
rela ted t o the difference in parameter space discussed in 
Sec. |3.I| Due to their exponential change in luminosity, 
the quasars are not typically found near their peak lu¬ 
minosities. Thus, even though by construction the peak 
quasar luminosity as a function of halo mass is the same 
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Figure 8. The halo mass function for halos hosting quasars versus 
the total halo population for certain models as defined in Table 
In the case of the lightbulb model, there is a roughly constant value 
of ~1% of halos hosting quasars as a function of halo mass. For 
the exponential, there is a constant active fraction at the high- 
mass end, but then the fraction falls off below the median halo 
mass. This is again related to the fact that there are more halos 
included at low mass in order to reproduce the average clustering 
signal. 


for the two models, the effective luminosity for a given 
mass is reduced in the exponential case due to the light 
curve evolution. In other words, quasars of the same lu¬ 
minosity in the two different models are found in more 
massive hosts in the exponential case. This leads to a 
systematic shift in the preferred mass range for cluster¬ 
ing measurements. To sum up: the increased spread in 
halo host mass for the exponential model compared to 
the lightbulb is due to inclusion of highly biased hosts 
in the measurement being balanced out by lower-mass 
ones, and the systematic shift toward higher mass is due 
to the effective increase in the mass-to-luminosity ratio 
related to evolution of quasar luminosity. 

4.2. Mass function and duty cycle of halo hosts 

To observe the effect that different points in parameter 
space have on the halo host properties, the mass function 
of halos hosting an active quasar has been calculated for 
the fiducial redshift selection. (For the hirfn and low- 
redshift selections, see Appendix [^) Figure]^ shows the 
total halo mass function as well as the mass function of 
halos hosting quasars within the luminosity range — 25 > 
Mi > —27. From this analysis, the duty cycle of halo 
hosts can be extracted, i.e., the fraction of active halos 
divi ded by the total number of halos. As discussed in 
Sec. |2.2[ in the lightbulb model the duty cycle can be 
directly related to the quasar lifetime at that luminosity. 
However, here the duty cycle is defined simply as the 
active fraction of halos. 

Figure [^compares the case of the lightbulb and the 
exponential light curves. When the halo mass function 
of active halos is examined in the two different cases, it 
can be seen that the mass range of hosts spanned by an 
individual model is quite different. In the lightbulb case, 
there is a very small range in host mass compared to 
the exponential case. This difference can be explained 
in terms of which hosts are included in the clustering 
measurements. In the lightbulb case, since the Inminos- 


ity is constant as a function of quasar lifetime, the only 
evolution in the relationship between mass and luminos¬ 
ity comes from mass accretion, which makes up a small 
fraction of total halo mass over the time scales for which 
quasars are active. As such, with an essentially static 
relationship, there is a very strong correlation of mass to 
light. For a specific model, there is only about a factor 
of 2 in halo mass included for the quasars in the selected 
magnitude range. When looking at the duty cycle of 
quasar hosts, one can see that the fraction is typically 
0.5-1%, with little evolution with mass within a model. 

Conversely, in the exponential model, there is evolu¬ 
tion for individual halos in the mass-to-light relationship. 
Most importantly, this implies that massive halos will be 
included when selecting quasars at a specific luminosity. 
Since they are more highly clustered (and more biased), 
smaller, less biased halos must also be included in or¬ 
der to create an average bias consistent with the BOSS 
measurements. This has the effect of extending the mass 
range of halos included in the mass function. Note that 
within a single model, there is a much larger span in halo 
mass: in some cases, the span is more than an order of 
magnitude in halo mass. Additionally, the duty cycle is 
comparable in magnitude to the lightbulb case, though 
slightly smaller: the ratio of active halos to total halos 
ranges from 0.05-1%. There also seems to be a trend in 
the evolution of the duty cycle: there is a central “typi¬ 
cal mass” for a given model, and the duty cycle de creases 
in both dire ctions. A similar trend was found by [White] 


etall(2012). 

JNote that one result of this measurement is the fact 
that the mass range of host halos is significantly more 
extended in the exponential case than the lightbulb case. 
Thus, one way to break the degeneracy between the light¬ 
bulb and exponential models would be to measure the 
mass range of underlying host halos, perhaps through us¬ 
ing gravitational leasing t o independently find the mass 
of the dark matter halo (Courbin et al. 2012|). If the 
range of masses for quasar hosts is extended, then there 
would be observational evidence favoring an exponential 
model (or a model with evolution in the quasar light 
curve) as opposed to the lightbulb model. 


5. PREDICTIONS FOR HELIUM REIONIZATION 


One very important prediction that we can make from 
our quasar catalog is the redshift of helium reionization. 
In order to understand in detail the implications for he¬ 
lium reionization, we need to run full hydrodynamic plus 
radiative transfer numerical simulations. However, we 
can perform a semi-analytic calculation in order to find 
a rough estimate of the redshift of reionization by com¬ 
puting the fraction Qi of the universe’s volume that has 
been reionized (also called the volume-filling fraction), 
where Qj = 1 represents a totally reionized universe (e.^., 
Madau et al.||1999} Furlanetto fc Oh|2008 ): 


dQi f d(l) - _ 

—— = / dL -—— - CaAUeQi, 
dt J r/He dL 


( 22 ) 


where nue is the number density of neutral helium, Ug 
is the number density of electrons, is the production 
rate of ionizing photons for an individual quasar, aA(T) 
is the recombination coefficient, and C = (r/g) / (r/g) is 
the clumping factor of the ionized IGM. The minimum 
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-1/912 normalization: 



Figure 9. The volume-filling fraction Qi of doubly ionized helium defined in Eqn. \22\ . In each plot, we show the fiducial values we 
have for Q, as a function of redshift, which has the parameters a = 1.7, (5 = 3, and normalizing the luminosity at 912 A following |Lusso| 
|et al.| (|2015ll. This leads to a redshift of reionization of z ~ 2.5, comparable to the redshift of z ~ 2.7 suggested by observation of the 
neiium Lyman-a forest. We show the change in Qi as a f unct ion of varying these parameters. Top left: we compare t he difference between 
using the composite QLF of SDSS-I-COSMOS (see Sec. |2.5| for more details) and the one in |Hopkins et al.| ||2007||. The shaded region 
reflects differences in ionization level due to jointly varying the parameters over the ranges specifie d m Table^ Top right: we change 
the UV SED of the quasar, which affects the normalization at 912 A. In addition to the SEDs from [Lusso et ah ( |2015^ and HRH07, we 
show the radio-quiet template from Shang et al. ( 201l| . Bottom left: we allow the EUV SED spectral index for A < 912 A to vary from 
1.4 < o < 2.0. Bottom right: we vary the clumping factor of the IGM, from 1 < (5 < 5. See the text for additional details. 


luminosity of the integral decreases as redshi ft decreases, IGM has been ionized, and all of the helium is singly 


et al. 1120121 IShen & Kelly||2012[ |Cen & Safarzadeh||2U15| vidual quasar N^. the SED of 

Lusso et al. 

(2015 

) is used 

Sijacki et al. 2015). Ihe clumping factor measures the ^.o convert the specific luminr 

sity at 2500 A to that at 


being averaged (or resolution in the case of simulations). 
Note that these calculations assume a primordial helium 
mass fraction of Yue = 0.24. Following the arguments 


in the appendix of Kaurov & Gnedin (2014), we choose 
the case A recombination coetticient, which assumes that 
photons emitted from recombination are not reabsorbed 
by a neutral atom, increasing the recombination rate. 

It is assumed that initially, all of the hydrogen in the 

Although the arguments presented in the cited work are in the 
context of hydrogen reionization, the same arguments can be ap¬ 
plied equally well to helium reionization. Essentially, the authors 


912 A. It is then assumed that quasars have an SED that 
follows a power law L^{u) oc v~°‘ for values of A < 912 
A. The fiducial value chosen is a = 1.7, also based on ob- 
servations of the r est-frame UV spectra of quasars from 
Lusso et al. (2015). This calculation includes all photons 
with frequencies m the range 54.4 eV <hv<\ keV. Pho- 

argue that the photons redshift out of resonance with the ther¬ 
mally broadened spectral line before they encounter the edge of 
the ionized region or a Lyman-limit system. Although the ioniza¬ 
tion fraction of helium might be slightly lower inside an “ionized 
region” than a comparable hydrogen one, the difference is not sig¬ 
nificant enough to change the overall conclusion. 
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tons above this energy have a mean free path of helium 
ionization comparable to the Hubble distance. 

Although the two different quasar light curves explored 
above have different individual properties, both are con¬ 
strained by the global properties fixed by the QLF. We 
find that if instead of the statistical calculation outlined 
above, we use the number of ionizing photons computed 
directly from the quasar catalogs, the result differs only 
by a few percent. The refore, it is much more straightfor¬ 
ward to use Eqn. (22). This approach also permits the 


use of other QLFs m the calculation, so it is possible to 
explore what effect this has on the results. 

Figure shows the ionization fra ction as a function 
of redshift computed from Eqn. (22). In the first panel. 


there is a comparison of the choice of QLF used in the 
calculation. Included are the QLF used in the main body 
of this work, the composite QLF composed of the ones 
from RI3, MI2, and MI3 (the SDSS-I-C OSMOS) as ex¬ 
plaine d in Sec. |2.5[ and the QLF from [Hopkins et al.| 
(2007) (hereafter referred to as HRH07). All other cal¬ 
culations presented use the composite QLF, but then 
change various other parameters. Note that for the pre¬ 
diction of reionization time using HRH07, both the QLF 
and the SED are different from the fiducial comparison 
case. In the figure, the shaded region shows the range 
of predicted values for the volume-filling fraction Qi at a 
given redshift z by jointly varying the parameters of the 
QLE over the range specified by Table Interestingly, 
the late-time ionization level is less sensitive to the varia¬ 
tion in parameters at early redshift, due to the interplay 
betwe en the source and recombination terms present in 
Eqn. (22 1 . At redshifts z < 3.5, the source term becomes 


the same for all histories, since the QLF transitions to 
that of Ross et al. Further, the recombination rate is 
proportional to the ionized fraction, so histories that had 
higher ionization levels at z > 3.5 will have higher levels 
of recombination. Since the recombination time is much 
shorter than the total timescale of the reionization cal¬ 
culation, all histories converge on a similar redshift of 
total reionization {Qi = I). Nevertheless, the variation 
in ionization fraction at early times can have important 
implications on the topology of ionized regions and the 
thermal history of the IGM, so such differences may in 
principle be detectable. 

In the second panel of the plot, the specihc luminos¬ 
ity of individual objects at 912 A L 912 is varied. One 
way to achieve this variation is the change the UV SED 
template used for quasars. Once the specihc luminosity 
7^2500 is calculated from the observed magnitude accord¬ 
ing to Eqn. (0, the quasar SED can be used to hnd L 912 . 
In the hducim app roach, we use the SED template from 
Lusso et al. (2015), which assumes a UV spectral index 
of a = 0.61 for 2500 A > A > 912 A. An alternative 
choice for an SED is one from 
provides a composite quasar S 


Shang et al. |20II), which 
iU template by combining 


Shang et al. (2011 1 divide the sample into 


radio-loud and radio-quiet quasars. However, radio-quiet 
quasars co mpose ^ 90% of h igh-redshift quasars found in 
the SDSS ( Shen et ar]|2009 ). Thus, we only include the 
results of the calculation using the radio-quiet template. 
This template provides the relative specihc luminosity 
at each frequency, and so can be used to convert T 2500 


to Lgi 2 . The effective spectral index for this wavelength 
range for the radio-quiet quasar template is a = 0.867. 
In addition, we show the impact o f usi ng the SED from 
HRH07 (with the QLE from Sec. 2.5). Note that the 
SED from HRH07 is outdated, and used only as a point 


of comparison. More r ecent studies {e.g., Stevans et al. 
20I4| Lusso et al.||20I5 ) are largely inconsistent with the 


SED of HRH07, and so it is presented here merely to em¬ 
phasize the importance that using the proper SED has on 
helium reionization. Given this same specihc lumi nosity 
^ 2500 , the p redicted value of L 912 from the SED ofjLussoj 


et al. (2015) is higher than that of HRH07 by about a fac- 
tor of 1.7, leading to the earlier reionization time. The 
second panel of the plot includes these to demonstrate 
the difference from using different quasar templates. 

In the third panel of the plot, the spectral indices are 
varied, ran ging from 1.4 < a < 2.0. Recent measure¬ 
ments from Lusso et al. (2015) suggest that at high red¬ 
shift and bright magnitudes, the spectral index has a 
value of a = 1.7 ± 0.6. This is slightly softer than the 
average value of a = 1.6 from Telfer et al. (2002). In 
order to explore some of the implications of changing the 
spectral index, we vary its value as indicated. 

The final panel explores a range of clumping values, 
from I < C < 5. The precise value for the clump¬ 
ing factor for helium reionization is very uncertain, as 
most studies on the clumping factor are related to hy¬ 
drogen reionization (see, e.^., 

Kaurov & Gnedin 2014). In 


Raicevic fc Theuns 2011 
I'urlanetto &: Uh ~(|ll008 ) 


the authors explored clumping factors of 0 < G < 3. 
More recent results from numerical simulations were cal¬ 


culated by Jeeson-Daniel et al. (2014), who found that 


the clumping factor of helium ranges from 3 < G < 8 for 
the redshift range of interest, depending on the ionization 
level of the helium gas. 

In Figure each panel shows the hducial evolution 
of Qi, which is char acterized by t he val ues of a = 1.7, 
G = 3, the SED of Lusso et al. ( 20I5|), and the com¬ 
posite SDSS-|-GOSM(JS QLI'. In this situation, the red¬ 
shift of reionization (ie., when Qi = I) is z ~ 2.5. 
This value is comparable to, though slightly later than, 
the redshift suggested by recent observations of z ~ 2.7 
(Dixon & Eurlanetto 2009 Worseck et al. 2011). How¬ 


ever, a smaller volume-averaged clumping factor G or 
a larger amplitude in either the measured QLF or the 
specific luminosity L 912 could give an earlier redshift of 
reionization. Specifically, assuming^ the hducial model, 
changing the clumping factor to G = 1.7 would give 
z ^ 2.7 as the redshift of reionization. It should be noted 
that this calculation is not wholly accurate for reioniza¬ 
tion, since it assumes a single clumping factor for the 
entire IGM, which is almost certainly not accurate for 
helium reionization, due to its very inhomogeneous na¬ 
ture. Furthermore, this calculation does not include sec¬ 


observations i n different frequency ranges to create a sin- 1979 
gle spectrum. 


ondary ionizations from energeti c electrons {e.g., Shull 
Furlanetto fc Stoever|2010), though these int erac¬ 


tions are likely uni mportant for helium reionization (Mc- 
Quinn et al.||2009 ). 


Whe n comparing to the results of [Furlanetto fc Oh 
(2008), we notice that the authors’ value for the redshift 


of reionization is significantly earlier than the one that we 
have found. This is largely due to a different QLF used, 
as well as a different method for calculating a quasar’s 
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EUV SED. The referenced paper uses the QLF from 
HRH07, and assumes an SED that gives more EUV radi¬ 
ation. This luminosity function has a significantly larger 
amplitude compared to the results from R13, up to an 
order of magnitude larger fo r low-luminosity quasars at 
high redshift. (See Fig. 16 of Ross et al. ]M^ ) Thus, ac¬ 
curate measurements and a proper understanding of the 
systematics of the high-redshift QLF, as well as the ac¬ 
companying quasar SED, are essential for a proper treat¬ 
ment of helium reionization. 


6 . CONCLUSION 

We have provided a technique for populating dark mat¬ 
ter halos with quasars that matches a QLF by construc¬ 
tion for various light cu rve models of quasars . By us¬ 
ing the triggering rate of Hopkins et al. (2006) with the 
technique of abundance rnatching, we are able to match 
the observed QLF of SDSS Data Release 9 (DR9) (R13), 
COSMOS (M12), and high-redshift SDSS data (M13). 
After applying this method to dark matter halo catalogs 
generated from V-body simulations, we have constrained 
a class of quasar models that reproduce the clustering 
amplitude measured from th e two-point auto-co rrelation 
function of the BOSS survey (|White et al.|2012|) at a red- 
shift of z = 2.39. The characteristic mass of the quasar 
hosts is 2.5 X 10^^ h~^MQ for the lightbulb model and 
2.3 X 10^^ h~^MQ for the exponential model. The effec¬ 
tive lifetime as defined in Eqn. (20) of quasars is teff = 59 


Myr for the lightbulb model of quasars and teff = 15 Myr 
for the symmetric exponential model. 

One of the limitations of this approach is that we have 
constrained the class of quasar models using a compar¬ 
atively narrow span in quasar luminosity. By matching 
the bias of quasars with a different magnitude range, we 
would have a different effective luminosity range for the 
bias calculation. This would lead to a different slope 
in the parameter Leff) which would allow us to break 
the degeneracy observed in Fig. Having the ability to 
break the sample down into different luminosity intervals 
would allow us to make tighter constraints on the class 
of allowed models. 

In future work, we plan to use the quasar models ex¬ 
plored here as sources of ionizing photons for studying 
helium reionization using simulations containing hydro¬ 
dynamics and radiative transfer. These types of sim¬ 
ulations will allow us to accurately capture important 
physical characteristics related to the IGM. Specifically, 
we are interested in capturing the thermal history of the 
IGM as it relates to observations. In upcoming simu¬ 
lations, we plan to compute the IGM equation of state 
and produce synthetic Lyman-a forest fluxes. This will 
allow us to tap into the wealth of observations available 
for the Lyman-a f orest, such as th ose currently available 
from BOSS {e.g., 
future surveys sue 


Lee et al. 20131, and from upcoming 
1 as DFSI - 
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APPENDIX 

A. FITTING THE PARAMETERS OF THE QLF 

In order to construct a QLF informed by the observa¬ 
tions at all redshifts relevant to helium reionization, we 
have combined the measurements of R13, M12, and M13. 
We will now briefly summarize the relevant findings of 
each paper. In all three results, the QLF is parameter¬ 
ized as a double-power law, according to Eqn. (II). RI3 


uses quasars identified from SDSS-HI DR9, and provides 
a LEDE model in which the base-10 logarithm of the 
QLF normalization, logj^Q 4>* j and the break magnitude 
M*, evolve lin early wit h redshift, as parameterized in 
Equations (12) and ([T^ . The parameters a and /3 are 
fixed as a function of redshift. Nominally, the LEDE fit 
is valid over the redshift range 2.2<z<3.5. M12 uses 
data from the COSMOS survey, and measures the four 
QLF parameters at z ~ 3.2 and z ~ 4. M13 uses quasars 
identified in SDSS data in Stripe 82 (S82), and reports 
the four QLF parameters at z ~ 5. For all three results, 
the parameters themselves and their associated la un¬ 
certainties are reported. In the M13 results, the authors 
actually provide three different fits to the observed re¬ 
sults. In their fiducial result, they fix the value of /3, 
and fit for the three parameters logj^Q (j>*, M*, and a. In 
a second set of parameters, the authors fix the value of 
a and find the best-fit values for the other three quan¬ 
tities. Finally, the authors fix M*, a, and P, and only 
fit for logigc^*. The best-fit values for the parameters 
change significantly in some cases between the different 
fits. More importantly, none of these fits seems to be 
ruled out conclusively by the data presented in M13, and 
so we incorporate all of the fits in our results. 

As explained in Sec. |2.5[ our goal is to combine the 
observational data froni different epochs. For redshifts 
z < 3.5, the parameters from R13 are used. At higher 
redshift, the parameters are assumed to vary linearly in 
redshif t. The eq uations for the parameters are given in 
Eqns. ( 15a|15d). The constant values are taken to be 
those ot Rid at z = 3.5, and the slope of the redshift 
evolution is allowed to take on a range of values. We will 
now discuss each of the four parameters in turn. 

For the parameter log^gc/)*, the fiducial value for the 
slope Cl is chosen to reproduce the average of the three 
reported values of M13 at z ^ 5. As discussed in M13, 
the fits from R13 extrapolated to z ^ 5 do not reproduce 
the overall normalization well, and predict too high a 
number density. Thus, a steeper value than that of RI3 
is necessary. The range of values for ci are chosen to 
bracket the range of best-fit values reported by M13. 


Trac, H., Cen, R., & Mansfield, P. 2015, ApJ, 813, 54 
Trac, H., & Pen, U.-L. 2004, New A, 9, 443 

Valageas, P., Clerc, N., Pacaud, F., & Pierre, M. 2011, A&;A, 536, 
A95 

Warren, S. J., Hewett, P. C., & Osmer, P. S. 1994, ApJ, 421, 412 
White, M., Pope, A., Carlson, J., et al. 2010, ApJ, 713, 383 
White, M., Myers, A. D., Ross, N. P., et al. 2012, MNRAS, 424, 
933 

Worseck, G., Prochaska, J. X., McQuinn, M., et al. 2011, ApJL, 
733 L24 

Wyithe, J. S. B., Loeb, A. 2003, ApJ, 595, 614 

Yu, Q., Lu, Y. 2004, ApJ, 602, 603 

Yu, Q., Sz Tremaine, S. 2002, MNRAS, 335, 965 

Zehavi, I., Zheng, Z., Weinberg, D. H., et al. 2005, ApJ, 630, 1 

Zheng, W., Meiksin, A., Pifko, K., et al. 2008, ApJ, 686, 195 






Figure Al. A plot of the evolution of the QLF parameters as a 
function of redshift: the base-ten logarithm of (p* (top left), the 
break magnitude M* (top right), the faint-end slope a (bottom 
left), and the steep-end slope f3 (bottom right). Best-fit values and 
associated Icr errors from R13, M12, and M13 are represented as 
the solid lines with shaded error regions, dark-gray triangles, and 
light-gray stars, respectively. For the M13 data, all three sets of 
parameters provided by the authors are plotted at 2 ~ 5, slightly 
offset for visual clarity. The dashed lines for 2: > 3.5 show the fidu¬ 
cial evolution of the QLF, and the dotted lines show the bracketing 
ranges of values explored. See the text in this appendix for further 
details. 


For the parameter M *, the fiducial value of the slope C 2 
is chosen to reproduce the average of the three reported 
values of MI3 at z ~ 5. The slope is allowed to take 
on a range of values that bracket the three reported val¬ 
ues of MI3. Also note that we have converted between 
magnitude systems using Mpz = 2) = M 1450 — 1.486, 
which assumes a spectral index a = 0.5. I f instead the 
value of a = 0.61 is used, as suggested by JLusso et al.| 
(|20I5|) and used in the calculations of Sec. then con- 
version is Mjjz = 2) = M^.^ n — 1.681. Further, if the 
SED from Shang et al. (2011) is used, the conversion is 
Mi{z = 2) = iWi 450 ~ 2T3y. The reason for the differ¬ 
ences is that the AT-corrections depend on the spectra l 
index of the SED (see Eqn. 3 offRichards et al. 2006). 
By extension, the QLF can be affected when cornbimhg 
different data sets. However, to be consistent with pre¬ 
vious works that have combined disparate data sets in 
this manner {e.g., RI3 and MIS), we use the conversion 
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given by assuming a = 0.5. 

For the parameter a, the fiducial value of the slope C 3 
is chosen to reproduce the average of the three reported 
values of Ml3 at z ~ 5. As with the other parameters 
discussed, a range of values is also explored which brack¬ 
ets all of the reported values of M13. Further, the value 
of a is bounded to lie where a > —2. For a < —2, the 
QLF does not converge for low-luminosity objects, and a 
cutoff luminosity must be specified below which quasars 
do not contribute significantly to helium reionization. To 
avoid defining such a cutoff luminosity, the value of a is 
bounded. As a practical matter, the ultimate goal of this 
project is to study helium reionization using full numer¬ 
ical simulations, where the minimum resolved halo mass 
will set the lower-limit of quasar luminosities. 

Finally, for the parameter /3, the fiducial value for the 
slope C 4 is chosen to reproduce the average of the values 
from M13 at z ^ 5. The range of slopes is chosen to 
bracket the values reported by M13. For the fiducial 
choice of slope, the value of /3 does not vary significantly 
with redshift. This range of values incorporates much 
of the parameter space constrained by M13, without the 
values of /3 becoming arbitrarily steep. However, the 
choice of /3 ultimately does not significantly affect the 
ionization level predicted by Eqn. (22). 

As a final note, the values of a ana at z ^ 3.2 from 
M12 are nominally inconsistent with the combined re¬ 
sults from R13. However, when looking at the results for 
individual redshift bins at z ~ 3.2 (e.(?.. Fig. 15 from 
R13), the uncertainties for the R13 values are signifi¬ 
cantly larger, and the results are largely consistent at la. 
The values of logj^Q (j)* and M* from M12 at z ^ 3.2 are 
consistent with the results from R13, and those at z ~ 4 
are consistent with the linear redshift evolution given by 
the re quire ment of matching the M13 data. Note that 
in Fig. |A1[ we do not plot the value of log;^Q at z ~ 4 
from M 12 , because the reported lower-bound of the er¬ 
ror bars is larger than the best-fit value, which must be 
positive. Despite this fact, the best-fit value is very close 
to the hducial linear evolution given here. 

Figure |A1| shows the measured parameters as a func¬ 
tion of redshift, as well as the assumed high-z evolution 
for each parameter. The solid lines and shaded regions 
show the best-ht parameters from R13, and the individ¬ 
ual points with error bars show the results from M12 and 
M13. The dashed lines show the fiducial choices for the 
parameters, which are chosen as outlined above. The 
dotted lines show the full range of parameters explored. 
The range of parameter combinations is applied to he¬ 
lium reionization in Figure]^ in the top-left panel. Note 
that, as discussed in Sec. this uncertainty primarily 
affects the early stages of reionization. Due to the re¬ 
combination term in the calculation of the volume-filling 
fraction and the fact that all reionization histories use 
the parameters of R13 at z < 3.5, the high-z values for 
the QLF do not ultimately affect the timing of reioniza¬ 
tion significantly; nevertheless, the different reionization 
scenarios can leave unique observable signatures on the 
IGM. 


alogs were also partitioned by redshift into a “high- 
redshift” and “low-redshift” sample in an analogous man¬ 
ner to the auxiliary BOSS samples. In the case of the 
BOSS results, the “fiducial” sample is actually the com¬ 
bination of the “high-redshift” and “low-redshift” sam¬ 
ples, so these two datasets are statistically independent 
of each other, but not the fiducial sample. For the pur¬ 
poses of comparing with the quasar catalogs, however, it 
is possible to compute ^(s) at distinct points in redshift, 
and compare with the BQSS results. The central red- 
shifts for the high-redshift and low-redshift samples are 
z = 2.51 and 2.28, respectively. Then an analysis similar 
to the above is performed, but at these additional red- 
shifts. This procedure yields further constraints on the 
bias as a function of reds hift in terms of the model pa¬ 
rameters to and 7 . Figure [BI] is similar to Figure]^ and 
shows how the selection of models varies as a function of 
redshift. In general, we find that the choice of parame¬ 
ters for our model to and 7 evolves slightly with redshift. 
In general, the BOSS measurements show an increase 
in bias with decreasing redshift. In order to accommo¬ 
date this increased bias, the model parameters must vary 
slightly. In general, the model favors quasars with in¬ 
creased lifetimes as redshift decreases. Despite this evo¬ 
lution with redshift, the relationship between log 4 Q(to) 
and 7 remains fairly linear, and it is still possible to pa¬ 
rameterize these models in terms of the characteristic 
lifeti me a nd luminosity factors tes and Leff as defined in 
Eqn. ( |20[) . 

Table 0 summarizes the changes in best-fit parame¬ 
ters as a function of redshift. Interestingly, these val¬ 
ues change somewhat: as structure continues to build, 
models with increasingly higher bias values are preferred. 
The fact that the best-fit values change demonstrates 
that the passive evolution of an increased clustering sig¬ 
nal within a given model is not sufficient; rather, this 
redshift evolution introduces additional constraints that 
we can use to select the most appropriate model. Never¬ 
theless, the results are consistent with no redshift evolu¬ 
tion. The results of White et al. (2012) also suggest that 
redshift evolution is minimal. Extending the clustering 
measurements to a larger redshift range could provide 
important constraints on the properties of quasar hosts. 


C. BIAS AS A FUNCTION OF LUMINOSITY 

We can also examine the dependence of bias as a func¬ 
tion of quasar luminosity. In the preceding analysis, we 
looked at the fiducial luminosity selection of the BOSS 
measurements for clustering, —25 > Mi > —27. In order 
to break the degeneracy in Fig. we explored the impli¬ 
cations of measuring the clustering of quasars with differ¬ 
ent luminosity cuts. We examined a high-luminosity cut 
Mi < —27, and a low-luminosity cut —23 > Mi > —25. 
Unfortunately, since the simulation volumes are only 
1 {h~^ Gpc)^, there are an insufficient number (~400) of 
high-luminosity objects to constrain the two-point corre¬ 
lation function. 

When fitting the functional form of the two-point cor¬ 
relation function, a power law is used: 


B. BIAS AS A FUNCTION OF REDSHIFT 

In addition to reproducing the “fiducial” sample from 

the BOSS results, the quasars from the constructed cat- n., , r , r , , 

T its the function are made tor cases where the exponent 
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Figure Bl. Parameter space evolution of fo and 7 from Eqn. as a function of redshift for the lightbulb model (left) and the exponential 
model (right). As redshift decreases, the space of preferred models shifts slightly toward those with higher intrinsic clustering. This is in 
addition to the passive evolution in clustering signal that each individual model experiences, which constrains the space of applied models 
somewhat. Nevertheless, the results are consistent with there being no redshift evolution. 



Figure Cl. The best-fit parameter sq for the two-point correla¬ 
tion function in the form g(s) = (s/sq)”^ as a function of power- 
law index 7 from Eqn. (|^ for the lightbulb and exponential cases 
using a fiducial (solid) and low-luminosity (dashed) luminosity se¬ 
lection. The gray shaded region shows the BOSS measurement for 
the fiducial luminosity cut. For the low-luminosity quasars, we see 
opposite trends for the two models. For the lightbulb, more neg¬ 
ative values of 7 mean that dimmer quasars have longer lifetimes, 
which combined with abundance matching implies they have more 
massive hosts. They therefore have larger values of sq compared to 
more positive values of 7 . In the exponential case, larger values of 
7 show more clustering because the bright quasars are longer lived, 
and are more likely to be included in the low-luminosity cuts while 
they are below their peak luminosity. Since they are abundance 
matched to more massive, highly clustered hosts, this leads to the 
behavior seen. See the text for further discussion. 


j3 is allowed t o vary, and o t hers with a fixed value of 
/3 = —2 as in White et al. (2012). In both cases, the 
clustering length sg increases tor larger values of the bias. 
To fit the best parameters, the parameters sq and j3 that 
minimized the value were found, where 

6 is defined as the difference between the average ^(s) 
and the functional form and C is the covariance matrix, 
calculated in the same way as in Sec. |3.2| These fits were 


made for the best-fit models defined in Eqn. (20) using 
the value s in Table lU 

Figure |C1| shows the value of the correlation length 
fits So for the fiducial luminosity cut —25 > Mi > —27 
(solid lines) and the low-luminosity cut —23 > Mi > —25 
(dashed lines) for the lightbulb and exponential models. 
The data are somewhat noisy, owing to the compara¬ 
tively large shot-noise error in the correlation function 
measurement. However, there does seem to be a trend 
emerging: in the lightbulb case, for more negative val¬ 
ues of 7 , the bias is larger, with the opposite trend for 
the exponential case. In the lig htbu lb case, this can be 
explained by noting, as in Sec. |3.2[ that in abundance 
matching longer lifetimes lead to a larger bias in the host 
halos. For negative values of 7 , less luminous quasars 
have longer lifetimes. Subsequently these quasars are 
being hosted in more massive halos. This means the clus¬ 
tering is stronger for large negative values of 7 , implying 
a larger value of sp¬ 
in the exponential case, the opposite trend is observed 
due to the presence of high-Lpeak interlopers. For pos¬ 
itive values of 7 , brighter quasars have longer lifetimes, 
and are more likely to be included in the low-luminosity 
selection. Since these hosts are abundance matched to 
occupy more massive, more clustered halo hosts, this 
leads to a stronger clustering signal, and a larger value 
of So- The evolution is not as strong as in the lightbulb 
case, however. In principle, the clustering measurement 
in different luminosity ranges could help break the de¬ 
generacy of best-fit models. 

Unfortunately, in practice this type of measurement 
might be difficult to actually make. The change in bias 
between the extreme values of 7 is not very significant, 
and the measur emen t is very noisy. The shaded gray 
region in Figure |CI| shows the current Icr bounds from 
the BOSS measurement, which has a larger spread than 
the variation in Sq as a function of 7 . Nevertheless, this 
ratio is a possible way to break the degeneracy between 
the different models. 












